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ABSTRACT 

In our current interpretation of the hierarchical structure of the universe it is well established that galaxies collide and merge 
with each other during their lifetime. If massive black holes (MBHs) reside in galactic centres, we expect them to form binaries 
in galactic nuclei surrounded by a circumbinary disc. If cooling is efficient enough, the gas in the disc will clump and trigger 
stellar formation in situ. In this first paper we address the evolution of the binary under the influence of the newly formed stars, 
which form individually and also clustered. We use SPH techniques to evolve the gas in the circumbinary disc and to study the 
phase of star formation. When the amount of gas in the disc is negligible, we further evolve the system with a high-accurate 
direct- summation A/-body code to follow the evolution of the stars, the innermost binary and tidal disruption events (TDEs). For 
this, we modify the direct A/-body code to (i) include treatment of TDEs and to (ii) include "gas cloud particles" that mimic 
the gas, so that the stellar clusters do not disolve when we follow their infall on to the MBHs. We find that the amount of stars 
disrupted by either inf ailing stellar clusters or individual stars is as large as 10~ 4 /yr per binary, higher than expected for typical 
galaxies. 



1. INTRODUCTION 

Super-massive black hole (MBHs) binaries are expected to 
form after major galaxy mergers. The main driving mecha- 
nism for the MBHs to sink to the centre is dynamical friction, 
where they will form a binary and start to shrink the semi- 
major axis on their way to the final merger. Slingshot of stars 
from the surrounding stellar environment help the binary to 
further decay by exchanging energy and angular mom entum, 
down to distances of about 1 pc ( Begelma n et al. 11980] ). How- 
ever, if the amount of stars to interact with is depleted, there is 
a risk of stalling, so that the MBHs would not coalesce within 
a Hu bble time. This is the so-cal led "last-par sec problem" 
(see Merritt & Milosavljevic 2005 for a review on the whole 
process and references therein). 

Key factors to surmount this last "snag" in the evolution are, 
among others, the fact that (i) in the case of binaries with a to- 
tal mass of < 1O 7 M , slingshot eject ions suffice to guarantee 
coalescence within a Hubble time (Milosavljevic & Merritt 
' 2003|) ; (ii) the role of gas may be crucial in the evolution of 



the binary, starting at larger scales. It might well be that in 
a merger of gas-rich galaxies, if MBHs are present, they will 
coalesce soon after the galaxies merge, in some 10 7 Myr, if 
the gas is distributed spherically. If, on the other hand, the 
gas is forming a nuclear disc, the galaxies need only to have 
1% of their to tal mass in ga s for t his to happen. ( [Escala et~aL| 
[2004} [2005] ). |Cuadra et al] ( [2009] ) found that such gas discs 
could indeed commonly help in the merger of SMBHs with 
masses in the range of our study, whilst this mechanism fails 
for masses larger than ~ 1O 7 M ; (iii) following with stellar 
dynamics, resonant relaxation creates a steady state current of 
stars which can be as large as ten times the non -coherent two- 
body relaxation (Hopman & Alexander 2006 ). This is a po- 
tential source of new stars populating the depleted loss-cone; 



1 Max Planck Institut fiir Gravitationsphysik (Albert-Einstein-Institut), 
D- 14476 Potsdam, Germany 

2 Departamento de Astronoirna y Astrofisica, Facultad de Fisica, Ponti- 
ficia Universidad Catolica de Chile, 782-0436 Santiago, Chile 

3 Max-Planck Institut fiir Astrophysik, Karl-Schwarzschild-Str. 1 , D- 
85741 Garching bei Munchen, Germany 



(iv) the work of Berczik et al. ([2006 ) shows that considering 
a non- spherically symmetric system the final parsec problem 
is largely solved (v) massive perturbers, such as giant molec- 
ular clouds or intermediate-mass black holes, can accelerate 
relaxation by orders of magnitude compared to two-body stel- 
lar re laxation, so that many new stars are supplied to interact 
with ( Perets & Alexander 2008 ); (vi) it has been observed that 
young, compact star clusters such as the Arches and Quin- 
tuplet systems reside near the Galactic centre. If these star 
clusters have masses larger than 1O 5 M , they can make their 
way down to the Galactic centre even if they start f rom a dis- 
tance as large as 60 pc within a few million years (McMillan 
& Por tegies Zwart 2003 ). The tidal stripping of these young 



stars could eventually provide the binary system with a new 
set of some « 10 5 stars; (vii) if intermediate-mass black holes 
(IMBHs), with masses ranging between 10 2_4 M. exist in the 
centre of clusters, it has been predicted that within the inner- 
most central 10 pc, we can expect to have some 50 IMBHs of 
masses 1 3 M Q , and still some of them at scales of a few mil- 
liparsecs ( Portegie s Zwart et al.|20 06 ). The interaction of one 
of these IMBHs with the binary of SMBHs would obviously 
accelerate the process of shrinkage. 

The studies just cited provide a number of mechanisms to 
make the binary shrink. We expect then that a typical binary 
will be able to reach sub-pc separations, especially in the case 
of relatively low-mass MBHs in gas-rich environme nts. In 
this study, we concentrate on such a case (see section [2J] for 
details), which is expected when the parent galaxies are gas- 
rich and large amounts of gas fall to the centre of the new 
system, together with the MBHs. At that location, the black 
holes get bound to each other, thus forming a binary, and are 



surrounded by a massive, parsec-scale gaseous disc (e.g., |Es- 
|cala et al.|20051|Mayer et al.|2007[|Dotti et al.|2007| ). 



Such gaseous discs are similar to proto-stellar discs: due to 
their high density compared to the central object tidal force, 
the discs will be locally unstable to self-gravity, meaning that 
perturbations in their density field will grow. However, if the 
gas is unable to cool efficiently, then the gas will not be able to 
contract and form clumps, and the density perturbations will 
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be sheared apart, creating a quasi-steady spiral pattern. Re- 
markably, the spiral pattern transports the angular momentum 
outward, making the disc behave as an accretion disc. On the 
contrary, if the gas is able to cool quickly enough, then the 
density perturbations grow and form clumps, which shrink 
and further accrete gas, breaking up the gaseous disc com- 
pletel y and turning it into stars - the so-called fragmentation 



(e.g.,|Gam mie 200 r||Rice et al.|2005[ [Nay akshin et al. 2007 
|Lodato|2007[|Alexander et al.|2008[|Paardekooper|2012| ). 

In either cooling regime, the situation where at the centre of 
the disc the central object is a binary will lead to a non-trivial 
interaction between them. On previous studies we have fo- 
cused on the inefficient-cooling regime, showing that torques 
between the gas and the binary will shrink the orbit of the 
latter, while the angular momentum is driven out through the 
disc ( [Cuadra et aL||2009l |Roedig et al^[20TTl [20T2| ). In this 
paper we present the first numerical study of the fast cooling 
regime in which the disc fragments into stars, and follow the 
dynamical evolution of the binary -stars system. We carry out 
our study in two stages (see also Khan et aL|2 012): first we 
model the fragmentation of the disc using smoothed particle 
hydrodynamics, and then we switch to our direct- summation 
A/-body models to both follow the long-term evolution of the 
system and to study the occurence of TDEs. 

The reason for this two-step approach is that we first need 
to model the gas hydrodynamics in order to follow the frag- 
mentation process of the gas, including the formation of stars 
and their growth via mergers and accretion of gas. In princi- 
ple, one could wait for the gas to disperse or be accreted, and 
simply continue the same integration to follow the dynamical 
evolution of the stars for long time-scales. However, the SPH 
code we are using is not designed to follow the collisional N- 
body dynamics of the system, therefore, it is necessary to use 
a different code that allow us to model the system of MBHs 
and stars in a meaningful way. 

2. FIRST STAGE IN THE EVOLUTION: GASEOUS DISC AND STAR 
FORMATION IN SITU 

2.1. Two MBHs and a circumbinary disc 

Following Cua dra et al.| ( |2009| ), we concentrate on a bi- 
nary with the following initial parameters: total mass Mbbh = 
M\ +M2, where M\ and M2 are the masses of the individual 
MBHs, and mass ratio M\ /M2 = 3, in a circular (Newtonian) 
orbit of separation a. The binary is surrounded by a corotat- 
ing gaseous disc with an initial mass M& = 0.2Mbbh and ra- 
dial range 2a-5a. The gas is modelled as an ideal gas with 
7 = 5/3, and radiative cooling is mimicked with a cooling 
time defined as t coo \(r) = f3/£l(f), where f3 is a free parame- 
ter that fixes the cooling rate, Q(r) = ^GM hhh /r 3 is the or- 
bital frequency around the binary, and r is the distance from 
the binary centre of mass. Since we are interested in the frag- 
mentation regime, in this paper we consider fast cooling rates, 
P < 5. The choice of a disc that fragments is realistic for self- 
gravitating discs t hat coo l therm ally, above a certain surface 
density threshold. |Levin| ( |2007| ) showed that, for the masses 
and distances we are interested in here, that threshold lies in 
the 10-100 g/cm 2 range. 

This model for the system dynamics is scale-free, mean- 
ing that it can be scaled up or down to different masses and 
lenghts. However, in order to introduce star formation and 
also to estimate the rate of tidal disruption events (TDEs), we 
need to choose physical units. With that aim we set the to- 
tal mass of the binary as Mbbh = 3.5 • 1O 6 M and we choose 



a = 0.04 pc. This would be a typical mass for binary black 
holes in t he range that could be detecte d by a LISA-like ex- 
periment ( [Amaro-Seoane et al.|2012a|b| ). The chosen separa- 
tion corresponds roughly to the value where we would expect 
binaries to spend the longest of their evolution in a simple 
model that considers binar y shrinking due to stellar sc atter- 
ing from a spherical cusp ( Milosavljevic & Merritt 2003) and 
torques from a non-fragmenting disc (see |Cuadra et alT[ 2009 , 

their eq. 12). 

While several studies (Amaro-Seoane & Freitag 2006; 
Amaro-Seoane et al. 2009, 20101 |Preto et al.||201l( |Khan| 
|et al.|201 1[|2012| ) have shown that stellar dynamical processes 
pump up the eccentricity of a binary MBH, in this case we are 
assuming the binary has reached the inner parsec in a gas-rich 
environment. In such a case, the dynamical friction of t he gas 
on the MBHs drives them to form a circular binary (e.g., Dotti 
|et al.| ( [2007] )). Thus we choose a circular orbit for the initial 
configuration. 

2.2. Implementation and treatment of the disc fragmentation 

To follow the process of circumbinary disc fragmentation, 
we use a modified version of th e smoothed particle hydrody- 
namic s (SPH) code (Gadget, ISpringel et al.|2001||Spr ingel 
2005), co mbining the numeric al methods of |Nayakshin et al. 
( 2007 ) and Cuadra et al. ( 2009). Here we only briefly describe 
the methods, and refer the interested reader to those papers for 
more details. We model the gaseous disc as an ensemble of 
initially « 2 x 10 6 particles of w 0.35 M each. The code cal- 
culates the gravitational and hydrodynamical interaction be- 
tween gas particles, plus the gravitational interaction between 
all particles, including the MBHs as well as the "proto-stars" 
and "stars" that form during the simulation (see below). We 
use a softening of 0.001a for the gas particles and of 0.01a for 
the proto-stars. The MBHs do not use softening, but a sink ra- 
dius within which gas particles are accreted. This radius had 
a value of 0.3a. 

As initial conditions, we take the initially-circular system 
modelled by |Cuadra et~aT| ( |2009| , at a time T « 500^. In 
this way we skip the transient initial evolution caused by the 
homogeneous initial conditions described in their work, and 
start from a steady-state configuration in which the circumbi- 
nary disc has developed spiral arms. Notice, however, that 
their simulations used j3 = 10, avoiding fragmentation. In our 
new simulations we set the value of f3 to either 1, 2, 3, or 
5. As a result, the disc now forms clumps, which grow in a 
runaway fashion. Treating this this with a pure SPH model 
is not feasible, as the growing densities require ever shorter 
time-steps. To circumvent this problem, we introduce sink 
particles to model the proto-stars that we expect would form 
in these large density regions. 

Proto-star particles are created when the gas density 
reaches 30 times the Roche tidal limit, M\ ) \^/(2irr > ). How 
many stars will form out of a gas density peak is a very com- 
plex question, whose solution is well outside the possibilities 
of our study. In our model we deal with this issue in an in- 
dividual particle basis, i.e., each gas particle is turned into 
one proto-star particle of the same mass. However, the newly 
formed proto-star particles can merge with each other, thus 
forming higher mass stars. The merger criterion is simply 
that their distance is smaller than 2f m R p , where R p is the size 

4 Notice that the choice of a = 0.04 pc is below the classical ~ 1 pc sepa- 
ration of the "final parsec problem", but for the range of masses considered 
in this work we deem it not a problem, as we summarized in the introduction. 
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of the proto-stars, which we typically take as 10 15 cm, and f m 
is a free parameter with fiducial value of unity that mimics 
the effect of gravitational focusing. The size parameter corre- 
sponds to ~ H/10 (where H is the disc scale-height), which is 
roughly the thickness of the gas arms we observe in the simu- 
lations. Thus, in a two-step process we are in principle allow- 
ing all the dense gas within the same overdensity to form one 
proto-star. However, we only allow the proto-stars to merge 
with each other as long as their masses do not exceed 30 M . 
Once they reach this mass we turn the proto-star particle into 
an actual star particle. The motivation for this threshold is 
twofold: numerically, we form an actual star out of ,> 100 
gas particles; physically, we avoid the rapid formation of ex- 
tremely massive stars. Stars can merge with proto-stars, but 
not among each other. 

Stars and proto-stars also grow by accreting their surround- 
ing gas. We use an Eddington-limited Bondi-Hoyle prescrip- 
tion to calculate their accretion rate, and then pick up at ran- 
dom enough particles from the (p roto-) star neighbour s that 



dom enough particles irom trie (p roto- ) star n eighbours that 
are merged with the sink particle ( Sprin gel et al.||2Q03] ). To 
calculate the Bondi-Hoyle and Eddington accretion rates, we 



-Hoyle and h,ddmgt< 
use the mass of the (proto) star, and a radius that is either the 
main sequence value corresponding to that mass (eq. 11 in 
Nayakshi n et al.|2007| ) for the star particles, or the fixed value 
R v for the proto-stars. This difference results in a much faster 
growth for proto-stars than for stars. 

The black holes also accrete the few gas particles that get 
too close to them. This procedure is done mostly to avoid the 
short time- steps that would be required to follow those gas 
particle orbits. Accretion on to the black holes is modelled 
simply with a sink radius - all gas particles entering the re- 
gion around 0.3a of either black hole are taken away from the 
simulation, with their mass and momentum being added to the 
corresponding MBH (Cuadr aet al.|2006| ). 

We have ran 6 different SPH simulations. Four of them use 
the fiducial values mentioned above, but differ on the strength 
of the cooling. We refer to these runs as betal, beta2, 
beta3a and beta5. Additionally, since we tend to form 
many very massive stars, we explore the effect of decreasing 
the numerical size of the proto-stars, hindering their growth. 
For (3 = 3 then we run two additional simulations, bet a 3b 
and beta3c, in both of which we use a smaller size for 
the protostars of 1436. 8 Rq instead of the fiducial value of 
14,368/?0. Run beta3c has however a larger gravitational 
focusing factor of f m = 10 instead of the fiducial f m = 1. For 
both extra simulations then there is a more severe (Edding- 
ton) limit on the accretion rate for the proto-stars than in the 
fiducial beta3a, while simulation bet a 3b has additionally 
a smaller likelihood of proto- stellar mergers. 

These choices in the conditions for gas cooling and for 
transforming gas particles into "stars" arguably capture a suf- 
ficiently broad number of potential fragmentation scenarios so 
as to envisage our analysis representative of a self-gravitating 
disc, within the limitations of the rather expensive numerical 
experiments. 



3. FRAGMENTING DISCS 

We run the SPH simulations of circumbinary discs for sev- 
eral hundred binary dynamical times. Due to the gas self- 
gravity, clumps grow in the disc. Given the short cooling 
times, these clumps contract, achieving the disc fragmenta- 
tion. In most simulations, after only « 2001^, the vast ma- 



jority ( £ 90%) of the gas is turned into stars, as expected. 
The system then reaches a quasi-steady state in which stars 
very slow ly accrete the tenous left-over gas (see Nayakshin 
|et al.|2007] ). The gas morphology at that stage for the differ- 
ent simulations is shown in Fig.[l] 

The fragmentation rate is set by the cooling time of the disc, 
thus discs with lower values of (3 will evolve faster. We can 
see this in figure[2] which shows the mass in stars as a function 
of time for all the simulations. The fourth column in Table Q] 
shows the number of stars formed in each simulation. Con- 
sidering only the variation of (3, it is clear that shorter cooling 
times result in la rger amounts of stars, as expected ( |Nayak-| 



shin et aO2007| ). As the total stellar mass is approximately 
constant, the typical stellar masses will be lower for shorter 
cooling times. 

It is interesting to note that the star formation process is not 
uniform. Instead, it happens preferentially in a few localised, 
relatively large regions, whose sizes are set by the spiral-arm 
overdensities. Even though we allow proto-stars to merge 
when they form close together, our numerical recipe avoids 
the formation of very large stars, which forces the formation 
of "stellar clusters" (see the left panel of Fig. [5| Some of 
these clusters feel a strong torque from the spiral arm and are 
driven towards the centre of the system, where the tidal force 
of the binary disperses them. This stellar distribution affects 
the long-term dynamics of the system and has interesting con- 
sequences for the production of tidal disruption events ( |4.1| ). 

In our tests with (3 = 3 and different stellar growth recipes 
we first notice that runs bet a3a and bet a3c are practically 
identical, and that run bet a 3b has the same curve of stel- 
lar mass growth. From this we conclude that in our simu- 
lations accretion is not important and that stellar growth is 
driven by mergers of sink particles. ^\ We also notice that the 
number of stars formed is about an order of magnitude higher 
in bet a 3b, which has 10 times smaller proto-stars than the 
fiducial run, and that the effect of having smaller proto-stars 
in the simulation is similar to having a shorter cooling time. 

To continue our study of the evolution of the MBHs and 
circumbinary disc system, we will take the masses, positions, 
and velocities of all sink particles and use them as input in 
direct- summation A/-body simulations. For simplicity, we 
take the snapshot at time T = 300 SIq 1 for all configurations, 
except for bet a5. Since in that run the evolution is slower, 
we use the snapshot at T = IOOO^q 1 , by which time 90% of 
the gas has turned into stars. 

4. THE ROLE OF STARS IN THE SHRINKING OF THE BINARY 

To analyse the dynamical evolution of the MBH binary em- 
bedded in the stellar system product of the stellar formation 
we use a direct-summation code, NBODY6. This is a very ex- 

5 For a movie of this simulation, visit the URL 
http : / /members . aei . mpg . de/amaro- seoane/ 
[frag menting - discsj 
The encoding of the movie is the free OGG Theora format and should 
stream automatically with a gecko-based browser (such as mozilla or 
firefox) or with chromium or opera. Otherwise please see e.g. http: 
|//en . wikipedia . o rg/wiki/Wikipedia : Media __help_ (Ogg) | 
for an explanation on how to play it. 

6 This is actually not surprising, as stars grow by mergers in roughly the 
dynamical time inside an overdensity, fdyn ~ (Gp)~ 1 / 2 , which corresponds to 
about hundred years for the density values required for the introduction of 
sink particles. On the other hand, the Bondi accretion rate for a solar mass 
sink, even for those very high densities, is only Meondi ~ lO -5 M0yr -1 in our 
models, so the time required to accrete a single SPH particle turns out to be 
hoc ~ 3 x 10 4 yr. 
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Fig. 1. — Gas density projected in the X-Y plane perpendicular to the angular momentum vector of the system for betal, beta2, beta3a, beta3b, 
beta3c and beta5, from the left to the right and from the top to the bottom. White dots represent the "sink" particles, i.e. the MBHs and the stars formed 
during the simulations. All snapshots are at T = 300 f^ 1 but for the last one, which was integrated up to T = 1000 Hq 1 , because in that run cooling is quite slow 
and the number of star s is still v ery low at earlier times (see Fig.[2j. Note that there is virtually no difference betwen beta3a and bet a 3 c. The figures were 
made with Splash (Price 2007). 
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Fig. 2. — Accumulated stellar mass formed in the disc in Mq for our 
fiducial case of a binary of 3.5 • 10 6 Mq . All simulations but for (3 = 5, which 
needs a bit longer, reach relatively fast the maximum of stellar mass and 
saturate with values below 1O 6 M0. 

pensive method because we integrate all gravitational forces 
for all formed stars at every time step, without making any a 
priori assumptions about the system. This code belongs to the 
family of dynamical codes for particle systems with relaxation 
processes of Sverre Aarseth. The code uses the improved 
Hermite integration scheme as described in (Aarseth 1999 



2003). Since these approaches integrate Newton's equations 
directly, all Newtonian gravitational effects are included nat- 
urally. More crucial for this subject is that it also incorporates 
both the KS regularisation and the chain regularisation, so 
that when stars are tightly bound or their separation becomes 
too smal l during a hyperbolic encounte r, the system is reg 



ularised ( [Kustaanheimo & S tiefel 1965). The advantages of 
this code as compared to the leap frog integrator of GADGET 
for our particular problem are obvious, namely the high accu- 
racy in the energy conservation, since we are interested in the 
correct evolution of the inner binary of MBHs as well as in 
potential TDEs. For this aim, as we describe later, we modi- 
fied the standard version of NBODY6. 

For each simulation, the initial masses, coordinates and ve- 
locities for the stars and MBHs are taken from the Gadget 
data at the times shown in table [T] At that moment, the gas 
mass - stellar mass ratio is very low (see table [TJ column 
^gas/M*)- The gravitational effect of gas is almost negligible 
and we do not include it in the simulations. Despite our limit 
to the growth of "proto-star particles" in the SPH simulations 
(see section [3]), some "star particles" did manage to achieve 
very large masses. We deem those unphysical, so in the initial 
conditions for our TV-body runs we replace stars with masses 
m abo ve 120M Q with a cluster following a Plummer distribu- 
tion ( |Plummer|1911| ) consisting of equal mass stars with total 
mass m and radius 



/ m \ 1 / 3 
R= V3M bb J r ' 



(1) 



with r the distance to the centre-of-mass of the binary. The 
last equation corresponds to the Roche lobe of the massive 
star with respect to the MBH binary with mass Mbbh- 

In our A/-body simulations, table [T] we exclude stars which 
are at a distance r > 100a, where a is the semi-major axis of 




Fig. 3. — Evolution of the density profile in one A/-body simulation 
betal at different times in the evolution. The dashed line corresponds ap- 
proximately to the position of the inner gap in the SPH simulation. 



Model 


SPH time 


M g as/M* 


NSPH 




Nsplit 


betal 


300 


3% 


2536 


1895 


4469 


beta2 


300 


7% 


1429 


1141 


2768 


beta3a 


300 


9% 


699 


585 


1924 


beta3b 


300 


9% 


5487 


4486 


5193 


beta5 


1000 


10% 


167 


144 


1146 


beta3b95 


95 




5540 


5540 


5540 



TABLE 1 

Initial data for the NBODY6 runs. Notice that we do not 

INTEGRATE RUN beta3c USING THE Af-BODY TECHNIQUE, BECAUSE IT 
TURNED OUT TO BE IDENTICAL TO beta3a. SPH TIME IS THE MOMENT 
AT WHICH WE STOP THE GADGET SIMULATION, IN UNITS OF Q^ 1 , 
Mgas/M* IS THE RATIO BETWEEN GAS AND STELLAR MASS AT THAT 
MOMENT, A^sph IS THE NUMBER OF STARS THAT HAVE BEEN FORMED AT 
THAT MOMENT IN THE GADGET SIMULATION, Af NB IS THE NUMBER OF 
STARS WITHIN A DISTANCE r < 100<2 FROM THE CENTRE OF MASS OF 
THE BINARY AND N split IS THE NUMBER THAT WE GET AFTER SPLITTING 

all very massive stars into sub -clusters, as explained in 
section[4] The reason why the last model has more stars 

THAN beta3b AT T = 300 IS BECAUSE IT CORRESPONDS TO A 

previous moment in the evolution and, as we explained 
above, protostars are allowed to merge with each other. 
This last case is a special one, and we ran a dedicated 

SIMULATION FOR IT. SEE SECTIQn [4~T1 ALSO, WE NOTE THAT WHILE 

the gas was originally distributed in a rather narrow radial 
range (2a-5a), we end up with stars even at distances > 100<3. 
This is due to N-body scattering, as many star particles are 
formed in relatively crowded regions of the disc. 



the MBH binary. We assume those stars would have only a 
negligible effect on the binary evolution. They correspond to 
about a quarter of all stars in each simulation. As shown in fig- 
ure this cut in the cluster did not affect its global structure, 
its density profile remains roughly constant at large radii. The 
figure also shows that the region inside a few times the binary 
semi-major axis gets depleted quickly by sling-shot interac- 
tions, as expected. 

In figure H] we see the evolution of a and e for all cases, 
integrated with NBODY6 with the results of the SPH simula- 
tions as input parameter. After some 10,000 orbits the binaries 
reach a stagnation point from which the decay becomes much 
slower. The decay rates (l/a)(Aa/ At) averaged over the time 
period from 0.5 Myr to 3 Myr are: betal : 7.2 x 10~ 9 yr _1 , 
beta2 : 4.0 x lO^yr -1 , beta3a : 8.0 x lO^yr -1 , and 
beta5 : 4.0 x 10~ 9 yr _1 (although for this case we start at 
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1000 ^q 1 , which means actually from 0.56 to 3.06 Myr). 
In the first 0.1 Myr of the evolution, the significant drop in 
semi-major axis corresponds to decay rates of 6.7 x 10~ 7 yr _1 
for betal, 5.4 x lO^yr" 1 for beta2, 5.4 x lO^yr" 1 for 
beta3a and 3.2 x 10~ 7 yr _1 for bet a5. 

The early dynamical evolution (first few hundred ^q 1 ) is 
dominated by close encounters between the MBH binary and 
stars on radial orbits (i.e. in the loss cone of the binary). This 
is naturally accompanied by a high rate of tidal disruptions 
(see figure [7]) and a strong change in orbital binding energy of 
the binary. In the following long-term evolution, the loss cone 
has been depleted and the binary is subject to the secular ef- 
fects of the disk as a non-spherical background potential. The 
effect of this type of mass distribution is a slow exchange of 
orbital en ergy but a rather efficient exchange of angular mo- 
mentum (Mer ritt & Vas iliev 2010), which is consistent with 
the significant increase in eccentricity that we observe in this 
phase compared to the very slow decay rates in the semi-major 
axis. 

4.1. An inf ailing cluster of young stars 

In the SPH simulations modelling disc fragmentation we 
see large amounts of stars falling to the immediate vicinity 
of the MBH binary. In particular, in simulation bet a 3b we 
observe an inf ailing cluster^] at T = 95 Qq 1 (see figure^. 

Since this is quite interesting from the stellar dynamics 
point of view, we run a dedicated simulation for this partic- 
ular situation with the direct- summation code. Nonetheless, 
at this early stage in the evolution of the disc, there is a signif- 
icant mass in gas which has not yet transformed into stars. If 
we ran the simulation without taking into account the gas, the 
small stellar clusters would dissolve, as their potential wells 
would be abruptly much shallower and the stars could not be 
held together. We therefore have to include a prescription in 
the A/-body simulations for the role of the gas, since includ- 
ing the gas particles directly is well outside the scope of our 
work. 

In this dedicated A/-body simulation we model each dense 
region that contains a non-negligible amount of mass as one 
particle with a big softening length. For this, we define a 
sphere at every region of interest. We then look at the SPH 
gas particle distribution and group together all particles within 
this region, compute their total mass, center-of-mass position 
and velocity and create one "cloud particle" with these prop- 
erties (see Fig. [6]). In the subsequent A/^-body simulation these 
particles are integrated separately, which required a modifica- 
tion in the code. In all gravitational interactions, the gravita- 
tional potential of the cloud particle seen by a regular star is 
then computed as 



GM C 

r c + e 



(2) 



where M c and r c are the mass and distance to the cloud par- 
ticle and e denotes the softening length, taken to be of the 
order of the size of the corresponding sub-cluster. Although 
the concept of cloud particles is already implemented in the 
standard version of NBODY6, we modified it to integrate the 
cloud particles taking into account the gravitational potential 
of the other clouds, stars and MBHs in order to follow cor- 
rectly the orbits around the central binary of MBHs. 

7 As in the former footnote about the movie, from T = 95 onwards in the 
simulation. 



Simulation 


TDEs (yr 1 ) 


betal 


1.1 - 10- 4 


beta2 


1.4- 1(T 4 


beta3a 


6-10- 5 


beta3b 


9 • 1(T 5 


beta5 


2-l(T 5 



TABLE 2 

Tidal event rates for the simulations of tableQ] 



The effect is that the particles in the sub-clusters now feel 
an additional gravitational force corresponding to the cloud 
and thus stay within their respective group for a longer time, 
which allows us to study their inf all and to analyse TDEs. 
However, after one close encounter of a gas cloud with one of 
the MBHs, the cloud would suffer a stripping from the cluster 
and now float around as an unphysically big agglomeration of 
mass. This means that we can get only a meaningful result for 
the very first encounter of each sub-cluster with the binary. In 
this respect, when estimating the TDEs for the infalling clus- 
ter, we will be giving a lower limit, since we cannot simulate 
realistically further interactions of the cluster with the MBH 
binary. In the right panel of figure [6] we show the distribution 
of stars in the X-Y plane after the first interaction. 

5. TIDAL DISRUPTION EVENTS 

During the direct- summation A/-body runs any star entering 
the tidal radius Rj of one of the MBHs is considered to be 
tidally disrupted and its mass is added to the m ass of the hole. 
For a solar-type star, this radius is (see e.g. Amaro-Seoane 
2012 for a derivation and examples) 



= R* 



1/3 



(3) 



In the last expression M B h is the mass of one of the MBHs, 
R* the radius of the star and m its mass. In order to estimate 
the radius of a star given its mass, we adopt the simple rela- 



tion R* oc m° 6 ([Demircan & Kahraman 



1991 



Gorda & Svech- 



nikov |199"8] ) with the normalization that a solar mass star has 
solar radius. Using this in Eq. |3j we can compute the tidal 
radius in solar radii: 



# T ,0 = 1.29m° 6 



Mi 



BH 



1/3 



(4) 



where m is the mass of the star in solar masses and the pre- 
f actor comes from an empirical fit for high-mass stars. 

In figure [7] we show the accumulated stellar mass fraction 
in tidal disruptions for all simulations. Based on this figure 
and for a time interval of 0.1 Myr after the initialization of 
the simulations, we can convert the values in tidal disruption 
events, as shown in table [2] Notice that the rate is actually 
much higher at the beginning of the simulations, but that result 
is likely to be affected by our initialisation choices. 

DISCUSSION 

In this work we have presented the first realisations of frag- 
menting discs around a binary of two MBHs in SPH with star 
formation followed by direct- summation A/-body simulations 
of the resulting systems. We have evaluated different frag- 
mentation scenarios based on an approximation for the cool- 
ing rate of the gaseous discs and different prescriptions for the 
growth of protostars. 
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Time (fl^ 1 ) Time (Q^ 1 ) 
^00 10000 20000 30000 40000 50000 Q ^00 10000 20000 30000 40 000 500 00 




— beta5 



0.5 1.0 1.5 2.0 2.5 3.0 ' 0.5 1.0 1.5 2.0 2.5 3.0 

Time (Myr) Time (Myr) 

Fig. 4. — Left panel: Evolution of the binary semi-major axis for all cases in the A^-body simulations in units of the initial orbital period Qq 1 and in Myrs. 
Case bet a 3b was not integrated for more than 10,000 orbits because of numerical issues due to the high number of stars in that simulation. Right panel: Same 
for the eccentricity. 





Fig. 5. — Same as figure[T]for the simulation bet a 3b. At the SPH time 95 the closest stellar cluster hurls itself on to the binary and leads to an enhancement 
in the TDEs. We take the position of stars and gas particles from the left panel to start a dedicated direct- summation A/-body integration which we name after 
this instant of time, beta3b95. 



When the gas is almost completely depleted, we take the 
masses, positions and velocities of the newly formed stars and 
feed them to the direct- summation A/-body integrations with 
the proviso that if the protostar has a mass above 12OM , we 
convert it into an agglomeration of stars following a Plummer 
profile of radius the Roche radius of the protostar to avoid 
artificially-created very massive stars. 

We find that the rate of decay in our direct A/-body simula- 
tions is slower than the w 10~ 6 yr -1 found in the SPH simula- 



tions of Cuadr a et al.| ([2009), when scaling to the same masses 
and separations. 

We simulate with a dedicated direct- summation integration 
the particular case of a simulation in which a cluster of stars 
that forms during the SPH simulation falls on to the binary, 
the case beta3b95. For this, we modify NBODY6 to in- 
clude "gas cloud" particles that allow the inf ailing cluster to 
hold together in the dynamical simulation in which we cannot 
realistically simulate the gas. 

Infalling clusters such as this and the scattering of isolated 
stars lead to a significant number of TDEs. To make an accu- 



rate estimation, we made a second modification of NBODY6 
to implement stellar tidal disruptions, and we find that the 
event rates lie between 2 • 10~ 5 - ~ 10~ 4 per system per year, 
which lies on the high side of current (uncertain) estimates for 
the TDE rate in standard galaxies, which typically lie between 



lQ^-lO^yr 1 (Phinney|l989, Magorrian & Tremaine|l999 



Syer & Ulmerll 1999b, and lie well within the observed rates 



D. 



onley et al.|2002[|van Velzen & Farrar|2012| ). A particular 
interesting signature of these TDEs is the reverberation map- 
ping" response of the circumbinary disc to a burst of emission 
produced by the TDEs. The light from the burst excites the 
gas in the disc, producing emission lines. The time- variability 
of the spectra, the echo of the TDE, during the months af- 
ter the burst could in principle allow us to constrain the disc 
structure (Brem, Amaro-Seoane, Cuadra & Komossa; part II 
of this paper to be submitted). 

While our simulations cannot follow the evolution of the 
binary for much longer times, it is interesting to ask the ques- 
tion whether the semi-major axis of the binary reaches dis- 
tances that would lead it to coalesce within a Hubble time 
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Initial positions of stars and gas clouds 



Positions of stars and gas clouds after T = 1300 yr 



3 0.0 
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-0.4 



Pi 
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0.4 
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x (pc) x (pc) 

Fig. 6. — Left panel: Initial configuration in the x-y plane perpendicular to the angular momentum vector of the system for simulation beta3b95 of table[T] 
Stars are shown with black dots, gas clouds with red circles. The MBHs are depicted with green circles. Right panel: The same system after the cluster falls on 
to the binary, after ~ 1, 300 yrs. Note the enhanced number of stars in the vicinity of the binary. This translates in a larger number of TDEs. 

1 me ( " } feet of the stellar disc is of about SE tot = GqM^ hh /a, where q 

is the mass ratio of ste llar mass to BH mass. Following an 
argument similar to e.g. Quinlan ( 1996|); [Sesana et aL| ( |2007| ), 
if we compare this to the orbital energy at semi-major axis a, 
one finds that the relative change after ejecting all the stars 
is SE tot /E « q/v, where v = 0.2 is the symmetric mass ratio, 
well below what would be necessary to shrink the binary by 
one order of magnitude. We note that indeed ejecting half of 
the stellar mass only shrinks the binary semi-major axis by 
< 25%, as we see in the first 3 Myr of our A/-body simula- 
tions, in figure [8] While this is true for our specific scenario, 
we note that further episodes of gas inflow towards the cen- 
tre could potentially trigger more episodes of star formation 
in the disc, which would lead to star scattering and a further 
skrinkage. Moreover, while we have focused on the effect of 
stars formed in- situ on the binary, but the system will be sur- 
rounded by a stellar cusp that constitutes an additional source 
of shrinkage for the binary. The supply of stars that will in- 
teract with it can be enhanced by additional mechanisms in 
a more realistic picture than that of an isolated, spherically 
symmetric galactic nucleus, as we discussed in the introduc- 
tion. 




Fig. 7. — Accumulated stellar mass in tidal disruptions for the different 
simulations of table ^ Since we only track them in the A/-body simulations, 
the curves start at T = 300 Q^ 1 but for model bet a3b 9 5, which as explained 
before, had a dedicated run. Because we cannot simulate realistically more 
than one infall of the cluster, we stop it after the first periapsis passage. 

because of the emission of gravitational radi ation, measur- 
able in a L ISA-like detector such as eLISA ( Am aro-Seoane| 
|et al.|2012a| . For this, the binary has to shrink from an initial 
semi-major of a ~ 0.04 pc down to a w 0.003 pc. This corre- 
sponds to an increase of orbital binding energy of about one 
order of magnitude. The net change in binding energy after 
an interaction with one bound star of mass m* can be esti- 
mated as AE* = Gm+Mbbh/a. We start the direct- summation 
simulations with a ratio of stellar mass to MBH binary mass 
of & 10%, so that ab definitio the stellar mass that is formed 
is not enough for the binary to shrink down to the phase in 
which the evolution is dominated by gravitational radiation. 
Indeed, if we consider all stars in the disc to be ejected, we 
estimate in the limit of this low mass ratio that the total ef- 
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Fig. 8. — Same as the left panel of figure|4]but including the number of 
stars in each simulation, for the same colour but in dashed lines. By the end 
of our numerical treatment we have lost at least 50% in all cases. 
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